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synopsis 


AVAIEESH KOMAR KHARE 
M.Tech. (Mech.) 

Indian Institute of Technology Kanpur 

KUTETO - ELASTODYH /MIC ANALYSIS OP SLIIER - CRAM 
MECHANISM "WITH PLEZIELY ATTACHED SHIER 

The mechanism examined is the common slider - crack 
mechaiism in which the slider mass is connected with Ifae coupler 
by means of a linear spring. The crank and the coupler are elastic 
and their masses and inertias have also been considered. Equations 
for Kineto-Elasto dynamic Analysis of the Mechanism have been derived 
and an iterative scheme is suggested to solve them. Priction is 
assumed small and the effect of clearances in joints and tolerances 
on links have been neglected. Some example problems are worked out 
and the results are compared with the available ones. 



CHAPTER I 


II® ROEOCT IOF 


1.1 General Introduction 

During last fifteen years a lot of work has been done in 
the field of design of mechanisms. Many investigators have contri- 
buted towards kinematic synthesis of mechanisms and it has proved to 
be a powerful design tool for engineers. Almost all the techniques of 
kinematic synthesis of mechanisms, work well and give acceptable re- 
sults for low-speed mechanism. However, in common, they suffer with 
one major saortcoming. This is the ''rigidity" assumption which pre- 
vails throughout tire literature with few notable exceptions. 

In high - speed mechanisms the inertia forces are large 
enough to cause speed - fluctuations as well as appreciable elastic 
deformations in the links of the mechanism, which can not be over- 
looked. The effect of inherent elasticity of a mechaiism system m^r 
not be pronounced in low speed light weight mechanisms, but in high 
speed mechanisms and in particular in high precision mechanisms, the 
analysis may not be perfect unless this effect is considered. 3?or 
example a kinematically synthesized gripp er mechanism considering 
all of its links rigid, may miss the target if it is run at a speed 
more than one - half of its design speed. 

In certain machines with reciprocating or oscillating 
motions proper timing of the motion is a prerequisite to smooth 
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operations. At high speeds owing to the large inertia forces? the 
main drive undergoes fluctuation of speed and the link - members of 
the mechanism are elastically deformed. This causes a difference 
between the actual kinematic characteristics of the mechanism at 
certain predetermined positions of the main drive and that predicted 
from the point of view of rigid body assumptions. 

The deviation of the mechanism characteristics from the 
expected ones perhaps motivated the workers to search for some new 
analysis aid synthesis techniques in which the effect of inherent 
elasticity of mechanism system was taken into account. The analysis 
procedures of the mechanisms have been grouped in the following main 
categories, j~11 "j 

Static Analysis : Determination of Inertia forces? stresses, str ains 

and deflections in members and joints due to 
external and / or gravitational loading. 

Kinematic Analysis : Examination of the displacement , velocity ratios, 

acceleration ratios etc., of a mechanism with all 
of its members regarded as rigid. The reference variable is usually 
a position parameter, 

Tynsmic Analysis s Determination of the displacements, velocities 

accelerations etc. of a mechanism, including 
derivations of inertia forces of a mechanism made up of rigid 
members. The reference variable is time. 
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Elastic .Analysis : Examination of stresses and deflections of si 

elastic system due to static load in order to 
determine system flexibilities or stiffness. 

El asto dynamic Analysis : Examination of displacements, velocities, 
accelerations , stresses, strains etc., of amoving elastic 
mechanism. Inertia forces are calculated by assuming all the members 
rigid. 

Kineto-elastodynanic Analysis (KEDA) : Examination of the displace- 
ments, velocities, accelerations, stresses, strai n s etc*, of 
a moving elastic mechanism. Effects of elastic deformation upon the 
inertia forces are included in the analysis. 

An ideal aialysis procedure of high speed mechanisms should 
not only consider the fluctuation# in speed of the main drive owing 
to the inertia forces and the inherent elasticity of the mechanism, 
but it should also take into account the damping characteristics of 
the system, clearances m the joints, tolerances on the links of the 
mechanism etc. However, such an analysis will be very much complica- 
ted and nobody hitherto has analysed the mechanisms considering all 
these effects together. However, attempts has re been made to dynami- 
cally analyse and synthesize mechanism with assumptions of varied 
degrees* In some cases particular links were considered to be elas- 
tic and masses to be concentrated at the hinges and so on. The 
analysis and ^mthesis of slider - crank mechanism seemed to have, 
judging from the published papers car this topic, attracted least 
attention of all the planar mechanisms. 
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1.2 Literature Survey 

LUrmg past few years valuable contributions have been made 
in the field of dynamics of mechanisms. Perhaps, due to the inaccurate 
response of the kinematically synthesized mechanisms at high speeds, 
the literature has recognized the need for dynamic analysis and syn- 
thesis techniques. However, in a great majority of work the dynamic 
analysis has been performed considering rigid links. 

Early work in dynamics of mechanisms concerned themselves 
with deriving the velocities, accelerations etc. of mechanisms. 

Quinn jj51 , 32*j has pioneered in the energy distribution method and 
has set the stage for further work. lor determining dynamic charac- 
teristics of mechanisms, Hii^chhorn {"*16 j has given the ’’rate of charge 
of energy approach” in which the time rate of change of total energy 
of the mechanical system is equated to the power input. Thus, if the 
power input function is Mown, this method will yield the time res- 
ponse of the system. 

Elastic deflections of significant magnitude in high speed 
mechanisms due to their inherent elasticity have attracted the wor- 
kers and various techniques have been developed to account for this 
effect. The use of springs in the mechanisms hare also been made in 
order to alter the dynamic performance of the given mechanism in 
high speed motion. This may improve the response* of the system, 
reduce the shaking forces transmitted to the fraae or support or 
produce a performance rtiich would otherwise be unsatisfactory with 
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the rigid member mechanism. Tan Sickle and Goodman 




used springs 


to adjust the speed of mechanisms. Imam j^19”| examined the response 
of a four - bar function generator with a longitudinal spring as a 
part of its coupler. Davidson 5 j has examined a slider - crank 
mechanism with the slider mass connected by a spring. In jj^cTj Benedict 
states that "the addition of springs is one of the easiest ways to 
balance torque and obtain a smooth - running machine" . 


Some authors 2, 3, 19, 20, 25 


have dealt with "elastic 
complex" systems i.e. a mixed mechanism in which few members are 
elastic and the remaining are taken to be rigid. Lagrangian mechar- 
nics or Energy methods are used to derive the equation of motion. 
Generally, as the solution of these problems is tedious, only one 
member is taken to be elastic and that too with only one elastic deg- 
ree of freedom. 

Dubowsky and Freudensteinj^6~j have made notable contribu- 
tion to the study of backlash aid have modeled an elastic mechanical 
joint, which they hare called as "impact pair", and have determined 
the dynamic equations for the relative motion of such joints. 

Slider - crank mechanism has been the centre of study of 
several workers. Jasinski, Lee and SandorJ^20, 2^| have studied the 
stability of this mechanism having a flexible connecting rod. 

Keubauer, Cohen and Hall 29~J have made a notable contribution by 
examining the vibrations of the elastic connecting rod of the slider — 


crank mechanism 
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liniecki 26 j considers the fluctuation of input crank 

-j. 

speed due to the effect of inertia forces. He deals with a procedure 
for ■the kinematic redesign of the slider - crank mechanism so as to 
minimise the error between the actual and the desired slider veloci- 
ties at three predetermined positions of the crank. 


Davidson 5j has exaained a slider - crank mechanism connec- 
ted in series with a spring and mass. He has taken crank and coupler 
to be massless and rigid and the only inertia load is due to the mass 
of -the slider. He h as also assumed a constant driving speed and with 
these assumptions he has synthesized the mechanism as a function 
generator. 

Sheerwood^30 considers the dynamic synthesis of a slider- 
crank mechanism so as to minimise the error between the actual end 
the desired time taken by the slider to travel a prespecified dis- 
tance. Sheerwood aid Hockey j^34 , 3Tj have considered the redistri- 
bution of coupler mass of a slider - crmk mechanism from balancing 
point of view. 


The literature reviewed above has not provided a method for 
"Kineto - Elastodynamic .Analysis’' (KEDA) of mechanisms in vhich all 
effects are considered. Only recently some contributions hare been 
made towards KEDA of mechanisms. 


Erdman and Sendor have called KEDA to be "a frontier 1 ' 
in mechanism design. Sadler and Sandorjj5!0 used harmonic analysis 
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to analyse the kineto - elastoctynamic deflections of the extension 
of the coupler of a four - bar path generator linkage. Erdman, 

Imam, Sandor and Oakberg j"s, 9, 10, l5|have dealt with a general 
method for KEDA and Kineto - Blastodynamic Synthesis (KEDS) of 
mechanisms using Kineto - El asto dynamic Stretch Rotation Operator 
(KEDSRO ) . They have started with the kinematic analysis to get 
system forces and elastic deformations, derive input link accelerar- 
tion using law of conservation of energy and get velocities and 
accelerations of remaining links from what they call ’’Kineto - Elasto- 
dynamic family tree”, obtained by differentiation of KEDSRO. The 


entire procedure is iterated till a suitable convergence criterion 
is satisfied. Imam j^18j deals with a KEDA and KEDS technique and 
without using ’’KEDSRO” synthesizes a mechanism with minimum weight 


criterion. 


1.3 Aim Of The Present Work 


As stated earlier kineto - elastodynamic analysis of high- 
speed mechanisms is a recently developed field of study and deserves 
attention of mechanism designers. Particularly m high speed high 
precision mechanisms (e.g. gripper mechanisms of computers and prin- 
ting machines etc.) the elastic deformations are large enough to 
affect the mechanism response and hence the desired accuracy can 
only be achieved if the effect of the inherent elasticity of the 
mechanism is also considered during analysis* 
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In tlie present work kineto - elastodynamic analysis of the 
eccentric slider - crank mechanism with flexibly connected slider 
has been performed. All the links have been considered to be elas- 
tic. The masses and inertia effects of all the links have also been 
considered. However, the effect of the clearances and tolerances, 
damping characteristics of the mechanism and the friction effects 
at the junction of two links have been neglected. 

The common slider - crank mechanism is a particular case 
of the mechanism considered here vhen the stiffness of the spring 
connecting the slider is considered to be very hi^a (theoretically 
infinite ) . 



CHAPTER II 


F0EMCJ1ATI0N OP THE PROBIB! AND SOLUTION TECHNIQUE 

2.1 Formulation Of The Problem 

Figure 1 represents the mechanism under study. The stiff- 
ness k of the spring connecting the slider with the coupler is 
assumed to be linear. The crank aid coupler are assumed to be 
linearly elastic and their elastic and section modulii are taken to 
be uniform aid constant. The mass of the coupler is concentrated at 
its centre of gravity G- (which is taken at its mid - way for 
simplicity.) The gravitational forces, being much small in compari- 
sion to inertia forces, are neglected. The joints are assumed to be 
frictionless and slider - friction has also not been accounted. 

There is no play in joints. linear superposition of rigid - body and 
elastic deformations is assumed. 

The mechanism driver rotates at an average angular velocity 

XI. 'The main drive undergoes fluctuations of its speed due to 

inertia of crank, coupler and the slider mass and at the steady — 

state the input crank angular velocity © is a function of input 

angle © only. The distance X represents the actual location of 

c 

1 

coupler end point C while X represents location of point D. TZhen 
the mechanism is motionless point D on the slider coincides with 
point C on the coupler. For very large values of k (a rigid connec- 
tion between slider aid coupler) C and 3) would coincide for all 



Figure showing slider-crank mechanism 
with f!e>nkiy attachfd sliier 



Figure showing element coordinates (6 ) 
and element forces (f) 




Figure sbowirg system cocrlinates t d) 
ar.d s> c tem forces IP) 



Jf!c 


Figure showing free body diagram of the 
coupler under the action of element forces 
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phases of the motion. But in the mechanism under consideration the 

slider is flexibly attached and has mass m^. Therefore poinr D is 

displaced by an amount (x - JC c ) with respect to the point C. The 

displacement X would be the algebraic sum of ’’rigid" displacement 
c 

and elastic deformation of the coupler end C in x - direction. 

2*1.1 The Analysis Procedure s 

Since crank and coupler are having mass and inertia and 
are elastic, it is not possible to obtain displacements, velocities 
and accelerations of different links at various positions of the 
crank directly* The displacements would constitute of two compo- 
nents - the dynamic displacements and the elastic deflections. The 
elastic deflections cannot be obt sine d directly as they are func- 
tions of system forces. Bor getting ’true’ elastic deflections, 
’true' system forces should be known. Hit the 'true* inertia forces 
depend on ’true* accelerations which in turn depend upon 'true* 
elastic def lections . Thus, system forces aid elastic deflections 
are inter - dependent and hence a suitable iterative scheme is re- 
quired to evaluate elastic deflections. The following procedure of 
analysis is adopted. 

The equation of motion of the slider and the energy ba- 
lance equation (the equation obtained try equating rate of dhange 
of total energy to the power input to the system) constitute the 
equations of analysis. The ordinary differential equations thus 
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obtained are non - linear aid coupled. The values of velocities and 
accelerations of centre of gravity of the coupler and displacement 
of the coupler - end G (fig. l) are required in order to solve them 
to yield the velocities and accelerations of the crank and slider, 
and the slider displacement. Since the required quantities them- 
selves are functions of velocity and acceleration of the crank, 
these equations cannot be solved directly. An iterative scheme is 
adopted to solve them. 

Firstly, all the velocity and acceleration terms required 
to solve the above two equations are kinematically expressed in terms 
of velocity and acceleration of the crank. The displacement of the 
coupler - end G is aLso expressed kinematically. The equations are 
then solved together using fourth order Range - Kutta method, yield- 
ing dynamic velocities and accelerations of the crank and the slider, 
and the slider displacement. The system forces are then obtained. 
Elastic analysis via flexibility approach yields element deflec- 
tions. Then the velocities end accelerations of elastic deforma- 
tions may be obtained. Knowing accelerations of deformation the 
dynamic accelerations aremodified. The modified system forces, ele- 
ment deflections etc. may now be obtained and the process is iterated. 
This scheme is iterated to desired accuracy so "that system forces, 
elastic deformations aid crank accelerations are best modified. 

Nov/ the modified values of velocities aid accelerations 
of centre of gravity of the coupler & ad coupler - end displacement 
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be obtained. With these modified values Hie differential equa- 
tions of analysis are again solved to give a better approximation 
for the velocities and accelerations of the crank and the slider and 
the slider displacement. The entire procedure is now iterated till 
the mechanism characteristics obtained from two consecutive itera- 
tions match to the desired accuracy. 

2*1.2 Derivation of Equations of Analysis 

Hefering to Pig. 1 the differential equation of motion of 
the slider may be written as 


Or 


* * 

m X + k (X - X ) = 0 

1 c / 

' * 2 2 

X + w X = w X 
n n c 


( 2 . 1 ) 


where 



Prom kinematic considerations, one gets 


and 

where 


Cos 

X 

c 




1/2 


r = radius of crank 

1 = length of connecting rod 

a = amount of offset 

© = angle of the crank with x - axis 
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0 = angle between the coupler and x - axis 

X = "rigid” or kinematic displacement of coupler end C, 
c 

It has been shown |j26j that if the sum of •— and ■“ 

* 

is an all (e.g. less than 0.3) 


1 - ( — Sin 0 + y ) 2 


1/2 


1 / r ^ a \2 
1 - - ( T son e + j ) 


with an, error of nearly 0.08^. 


Therefore, X can be rewritten 


X 


r Cos © + 1 


-i {! s “ e+ if 


( 2 . 2 ) 


The working section of the torque - speed characteristics 
of a three phase a/c motor may be successfully approximated by a 

parabolic equation when referred to the main shaft p 26^J . 

.2 

T = — &2 ® 


(2.3) 


where 

T 

§ 


external torque overcoming fluctuations in speed 
referred to the main shaft. 

Coefficients of applied torque 
Crank velocity 


Applying the Energy Balance to the mechanism system yields 

d ( T.. E.) = T d 6 (2.4) 

mere, T.E. is the total energy of the system at any phase and 
includes the kinetic energy of crank, coupler and slider mass, the 
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potential energy stored, in the spring and the Strain Energy stored 
In crank and coupler owing to elastic deformations. But due to the 
presence of a spring, the strain energy stored in the crank and the 
coupler vail be much smaller in compamsion with the kinetic and 
potential energies and hence is not accounted. 

Therefore, 

T«- E. — K. E. + P. E. 

where, 

K. E. = kinetic energy of the mechanism 
P. E. = Potential energy stored in the spring 

Substituting in equation (2.4) results 

d (K.E. + P.E.) = T d 0 


(K.E. ) 




T = 


( 2 . 5 ) 


The K.E. of the mechanism is given by 


K.E. = 


I © + I *0 + m 2 + m 2 ^C- + m i 

O Gr 


(2.6) 


where 

I = Moment of inertia of crank about the axis of 
o 

rotation 

I = Moment of inertia of coupler diout its centre 

Gr 

of gravity 

jZ? = Angular velocity of rotation of the coupler 

X„ and Y„ = x and y - directional velocities of centre of 
G G 

gravity of the coupler 
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X » linear velocity of the slider 
m 1 = Mass of the slider 
Dig = Mass of the coupler 


Differentiating equation (2.6) with respect to 8 and setting 

* 2 

Y = 9 , one gets 


i <*•*•) - f 


i Q H0 + m 2 (z ! . X G+ T 0 y g ) 


— EL 
2 1 


t P ' I I! ■ 

X Y + 2Y X X 


( 2 . 7 ) 


Since 


i (X 2 ) = x' 2 § + 2Y x' X W 


and 


d9 


< 3 . t )) 2 \ 2 0 0 

r \ 0 ; = — and so on. 
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The primed quantities in the above equations represent 
derivatives with respect to 0 . 

The potential energy stored in the spring is given by 

P.E. = ~ % ( X - X Q ) 2 

Differentiating (2.8) with respect to 9, one gets 

(P.E.) = k (X - X c ) (X* - x ' 0 ) 

Substituting (2.7) and (2.9) in (2.5) and simplifying, one gets 

1 2 

, r I m X -r , „ 

Y I “g 2 + j= a., - Y (ag + m 1 X X ) - 

- (} S * V + -2 + \ V] 

- k ( x - x_) (X 1 - x' 5 


( 2 . 8 ) 


( 2 . 9 ) 


( 2 . 10 ) 
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Since 

X 


*2 n •• i 

e x + © z , 


equation (2.l) yields 

= Y < W n ^ 2 C- X )- (2.11) 

Prom (2.ll) aid (2.10), one can get 


» 

Y 


2_ 

I 


a. - aj Y - 


J T 


+ m 2 (i G ■£ + i V ) 1 


1 . 


- k (Z - X G ) (X* - X G ) + z' w^ (X- Z c ) 


( 2 . 12 ) 


Equations (2.1l) and (2.12) constitute the equations of analysis and 

are two non - linear coupled ordinary differential equations whose 

• « • * 

solution would give the values of X, X, ©, X and © . 


2.1*3 Elastic Deformation Of The Mechanism 


Any mechanism may be considered a stricture if its rigid - 
body - kinematic degrees of freedom are removed. At every phase the 
mechanism is "frozen 1 ' by removing its rigid - body - kinematic degrees 
of freedom and then its analysis as a structure yields the elastic 
deformations. 

A mechanism is composed of various combinations of elements 
each of which can be represented by a known structural model. The 
deflections of the entire mechanism system may be derived by per- 
forming an elastic analysis using flexibility approach. The mecha- 
nism will have system, or generalized external forces acting upon it 
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*1 a = i, 

3 i 


m, 


vhich will be represented by Hie column matrix 
viiiere m is the number of system forces. The number of the elastic 
degrees of freedom n of Hie system is the sum of elastic degrees 
of freedom of its elements, each degree of freedom being represented 
by an element coordinate. For the mechanism under consideration Hie 
number of elastic degrees of freedom n is taken to be 7 while 5 
generalized forces act on the mechanism. (Fig. 2 and Fig. 3) 


2. 1.3.1 Matrix Of Element Flexibilities 


The mechanism constituting of craik aid coupler only is to 
be analyzed for elastic deformation at each phase. It has got only 
one rigid - body - kinematic degree of freedom. Therefore, it may 
be converted into a structure by modeling the input link as a canti- 
lever or fixed - free beam. As the end C of the coupler is assumed 
to move along x - axis only, the coupler is modeled as a beam hinged 
at the end A and having double roller support at the end G as shown 
in Figure 2. 


*[*] 


r. 


The element flexibility matrix. F j of the crank is given 

L C i # 


F 

i cr 


r/A Q E 1 
0 
0 


0 


rV3E 1 I 1 
r/2E 1 I 1 


0 


r /2E 1 I 1 

r/E i h 


( 2 . 


'inhere 


A = Area of cross - section of the crank 
o 

= Modulus of elasticity of the material of crank 
1^ = Section - modulus of the craik. 


13) 



20 


The element flexibility matrix of the coupler cm be cons- 
tructed by using C astigliano ' s Theorem. 

Refering to FIG. 4 the strain energy per unit volume stored 
in the material of coupler is given by 


1 I f l/2 (f 4 Cos # + E c Sin 0) 2 

— I - — - dz 


T (f . Cos 0 + R_ Sm 0 - If. Cos 0 + f „ Sin 0 ) 

+ 4 C „ . 6 7 - dz 


V 2 (Py, Cos 0 - 1 . Sin 0 x2 


E 2 J 2 


2 A 

z dz 


f 4 z Sin 0 + R^ z Cos 0 + f g (z - -|) Sin 0 + f ? -Gz - *|)Cps-:-0 - f $ ) 


E 2 I 2 


( 2 . 14 ) 


IShere 


R c » the 


- Area of cross - section of the coupler 

E 2 = Modulus of elasticity of the coupler 

Ig = Section - modulus of the coupler 

R = Reaction at C 
G 

reaction at the point C, can be evaluated by taking moment of 


all the element forces about A. Then one gets 


= - 

c 


f 5 f 7 

Tosrr* f 4 


tan 0 + 


( 2 . 15 ) 
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Uow applying Castigliano' s Theorem, the differentiation of equation 
(2.14) with respect to f^ results 




a*. 


(f^ Cos 0+Rq Sin 0) dz 


. v 'o 

J. 


E 2 A* 


Cos 0 


r (f^ Cos 0 + Bq Sin 0 - f g Cos 0 + f ^ Sin .0) 

J e 2 X 

1/2 2 ^ 


dz 


Cos 0 


viiich results after simplification as 


1 Sec 2 0 tan 0 Sec 0 1 Sec 2 0 

4 4 Eg 5 B 2 A? 6 2 E 2 Aj 

Differentiation of equation (2.14) with respect to fg results 


(2.16) 


6, 


b U 


+ 0 


1/2 


^ 2 (f 4 Cos 0 + E c Sin 0) tan 0 


J, 


E 2 ^ 


dz 


(f^ Cos 0 + E q Sin 0 - f g Cos 0 + f^ Sin 0) tan 0 


E 2 A C 


dz 


f ^ 2 (E c Cos 0 - f 4 Sin 0) z 2 




E 2 J 2 


dz 


(“f^ z Sin 0 + z Cos 0 + f g (z - -|) Sin 0 + f ^ (z - ^)Cos 0 - f^) 


1/2 


E 2 ^2 


( f - 1 ) dz 


which, on simplification, can be written as 



22 


£ = ± Secjg tan 0 

°5 *4 e 2 5 


1 t an 0 

12 Eg Ig 1 Eg 


- £ , 


Sec 0 tap. 0 

2 


6 2 Aj, 


(2.17) 


Eiff e £e&tiat:ing equation (2.14) with respect to f fi , one gets 


£ , ^2. 


= -j 


f (f. Cos 0 + B. Sin 0 ) Sin 0 tai 0 


E 2 A} 


dz 


-i 

1/2 


1 (f 4 Cos 0 + 1^, Sin 0 - f g Cos 0 + f ? Sin 0) 


E„ 




(Cos 0 + ) dz 


^ (- f 4 z Sin 0 + z Cos 0 + f g (z - *|) Sin / + f ? (z - -|)Cos 0 - f 5 ) 


J, 

1/2 


®2 ^2 


( z . t ... ! , ) , cjjj 0 dz - 


1/2 


(e c Cos 0 - f^ Sin 0) Sin 0 g 


Eg Ig 


■£ — z dz , 


Or 


f _ _ -p 1 Sec 0 Sec 0 tan 0 

°6 “ *4 2 Eg *5 2 Eg ^ < 


1 (Sec 2 0 ± Cos 2 0 } 

4 E 2 Aj 


+ . £ su sljL I + f 

^Eglg i ' 


l 5 Sin 2 0 _ 1 tan 0 (l + Cos 2 0 } 

2 I 2 8 B 2 


96 E„ I n 8 E 0 _J 


( 2 . 18 ) 


Differentiating equation (2.14) with respect to f^, one gets 
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(f^ Cos 0 + Eg. Sin 0) Sin 0 


E 2 Aj 


^ (f^ Cos 0+3^, Sin 0 - f g Cos 0 + f ^ Sin 0) Sin 0 


(^j Cos 0 - f .Sin 0) 2 

rspr 2 2 Oos 0iz 


y (- z sin 0 + z Cos 0 + fg (z - Sin 0 + (z - 


■g)Cos 0 - f 5 ) 


*2 


Cos 0 dz 


Simplification of the above equation results in 


I 2 Sin 2 
6 96 B 2 I 2 


8E 2 V 


3 2 


4. -f . r. . ■arm , + -s — c 

T 4 E 2 48 E 2 X 2 


(2.19) 


The element deformations axe represented by an — dimensional 


column matri 




and are related to the element forces 


K) 


by the 


relation 


t - ki -? 


( 2 . 20 ) 


■where j""f^~j is n - dimensional column matrix representing the 
element forces and j~ PjJ is the element flexibility matrix (n x n) 
of the mechanism and can be constructed by using equations (2.13)* 


(2.16), (2.17), (2.18) and (2.19) as 
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E 


0 

0 

0 

0 

0 

0 


0 

0 

0 

0 

0 

0 

°2 

°3 

0 

0 

0 

0 


°4 

0 

0 

0 

0 

0 

0 

°5 

°6 

°7 

0 

0 

0 

°6 

C 8 

°9 

0 

0 

0 

°7 

C 9 

°10 

C 11 

0 

0 

0 

0 

°11 

C 12 


where 


C. = 


A o E 1 

r? 

3 E I 
0 1 1 


C„ = 


2E 1 X 1 


c, = 


cl = 


E 1 E 1 


1 Sec 0 

E 2 A} 


0, = 


tap. 0 Sec 0 

E 2^ 


C„ = - 


1 Sec g? 
2E 2 ^ 


C 8 = 


12 E 2 I 2 


tan 2 # 
1E 2 % 


( 2 . 21 ) 


_ Sec 0 tan 0 

9 "~2E 2 ^ 
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C - 1 (Sec 2 gf + Oos 2 ef) l 3 Sin 2 gf 

10 4 E 2 Aq 48 E 2 I 2 

n _ l 3 Sin 2 g f 1 Sin 2 ef 

11 " 96E 2 I 2 ~ SE,^ 

c - 1 Sin 2 0 , l 3 Cos 2 ef 

12 - 4^^ 48E 2 I 2 ( 2 . 22 ) 


2.1. 3.2 Eorce Transfer Matrix 


In order to transfer the system forces j P.l ] = 1, 

; D * 


• • • J 


m 


into element or internal forces j~f . j i = 1, . n an n x m force 
transfer matrix jjDj is derived by the methods of static analysis. 


The transformation relation is given by 


= ] 
_ 


D j P 

J 


(2.23) 


The force transfer matrix JjdJ also relates the system deformations 

represented by m - dimensional column matrix | d j and element 

J i 


deformation 


sn 


and the relation is given by 

~T = [d] 1 gT (2.24) 

Equations (2.20), (2.23) and (2.24) can be used to get the relation- 


ship between system deformations id. \ and system forces 

L il 

Thus substitution of (2.2o) and (2.23) in (2.24) results in 

T 


P. 

0 


or 


T - M N l »3 * 

■ 0 J r 


~T 


(2.25) 

(2.26) 


where 
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is ai xi matrix and is known as system flexibility matrix. 

!Ehe force transfer matrix {jdJ can be constructed by using 
the methods of static analysis. 

Refering to FIGURES 2, 3 end 4 representing force diagrams 
of the elements, one gets 


f„ = P 


1 


f c 


f, = P, 


= P, 


( 2 ; 28 ) 


« H Cos 0 - E Sin 0 
a a 

~ H Sin 0 + R Cos © 
a a 


(2.29) 


as 


The system forces cm be written xn terms of accelerations 


** 1 t 

I B = ■— I Y 
o 2 o 


h* 


m 2 X G 


(2.30) 


‘4 


k (Z - X 0 ) 


ss m« Y„ 
2 G 


The reactions H and R caa be obtained using free- body 
a a 

- diagram of the coupler RIG. 4. Thus, 


H 


f. - f, 
4 6 



f 6 f c f 7 

P a = 2 "k* 31 ^ + ~ ^ os ^ + ^4 ^ ai $ + ' 
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( 2 . 31 ) 


Substituting (2.3l) ixl (2.29) and simplifying, one gets 

■p _ _ -d Sin. Q . „ f Sin © ten 0 „ ^ 

f 1 “ P 2 1 Cos + P 3 ( 2 “ Cos ®) 

P- 

+ P4 (Cos 6 - tan jZf Sin ©) - -y- Sin © 


and 

f- 


Cos © 


2 1 Cos 0 ~ P 3 (Sin © + 


Cos © taa 


4 


JtV 

+ (Sin © + tan 0 Cos ©) + — Cos © 
from (2.28) and (2.32), one can obtain that 


( 2 . 32 ) 


D = 


where, 


0 

0 

1 

0 

0 

0 

0 


» n = 


D, = 


D 1 

D c 

0 

0 

1 

0 

0 


D 2 

D 6 

0 

0 

0 

1 

0 


D 3 

®7 

0 

1 

0 

0 

0 


D = - 


- Sin © 

1 cos 0 

s . tojLtgi J _ 0os e 

Cos © - ten ^ Sin © 

sin © 


D, 


8 


0 

0 

0 

0 

1 


(2.33) 


'2 



Cos 0 
1 Gos 0 


28 


D 5 - 


D 6 = - (Sin 8 + C ° B | tal 
D 7 = Sin 0 + ton 0 Cos 0 


D S " 


Gos 0 
2 


2. 1.3.3 System Flexibility Matrix 

Finally the system flexibility matrix, an 
mqy be derived by substituting equations (2* 21 ) and 


F 


where 




P 2 = C 3 *>5 


*3 " G 3»6 
*4 = °3 D 7 


'l 

*2 

*3 

*4 

*5 

2 

*6 

F 7 

*8 

P 9 

*3 

*7 

P 10 

F 

11 

F 

12 

P 4 

*8 

I.. 

11 

F 

13 

F 

14 

F 5 

P 9 

F 

12 

F 

14 

F 

15 


F r 


° 3 ^ 


"6 - °1 D 1 + C 2 D 3 + °8 


f ? = ° 1 d 1 d 2 ^ 2 d 5 d 6 + c 9 


(2.34) 

(m x m) matrix, 
(2.33) in (2.27). 


(2.35) 



29 * 


p. 


8 


G 1 I> 1 D 3 + G 2 I> 5 D 7 + Gg 


p. 


P 5= 
10 


p = 
11 


P a 
12 


P as 

13 


p = 

14 


P 


15 


°1 B 1 B 4 + C 2 B 5 S 8 
°1 + °2 D 6 + °10 
°1 B 2 D 3 + °2 B 6 ^7 + °7 
°1 B 2 B 4 + °2 B 6 B S + °11 
°1 B 3 + °2 B 7 + °5 
°1 B 3 “4 + °2 B 7 B 8 
°1 D 4 + °2 4 + °12 


(2.36) 


2. 1.3*4 System deformations end their velocities end accelerations 

Having derived the system flexibility matrix, the system 

or generalized deformations jjfTj j = 1, m may be obtained 

directly from equation (2.26) 


Eefcring to PIG. 3 one can conclude that 

e = d, 

e 1 


= - 2 


X, 


'G. 


X, 


« <1/ 


G 


( 2 . 37 ) 


where suffix e denotes the elastic deformations (e.g. 2L is the 
elastic deformation at point G of the coupler in x - direction)^ 
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The velocities and accelerations of defomaticn nay then, 
he obtained by using finite difference technique. Thus, one nay 
derive 


©. = 


0 - 


© & 

e 

t ♦ 

9 


X„ = 


X G 9 
e 


y g - 


y g 9 

e 


X„ = 


x e 

e 


( 2 . 38 ) 


and 


e 

• • 


e 


» 


n *2 


e ~ e e * e e « e e + e i 


jzf - 




y g - 


x„ = 


,1 ** It 

0 e e + 0 e -i 

I ** II 

Xq & + X G Y 

e e 

i ** it 

Y s 9 + Y g Y 

e e 

I *• N 

\ e * V 


(2.39) 


2.2 Solution Technique 


After having mathematically formulated the problem, the 
next task is to devise some iterative technique to solve the equar- 
tions. As discussed earlier the equations of analysis cannot be 
solved directly. 

Rofering to equations (2* 1 1 ) and (2.12), one may con- 
clude that this set of two non - linear coupled ordinary differential 
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• • 

equations may be solved directly to yield the values of Z, Z, 9, 

•• *• 

Xf © ote* provided the displacement, velocity and acceleration 

* * * 

terns (o.g. Zq, Z^, Z^ctc.) appearing an the right hand side of 
the equations are known. But these displacement, "velocity aid 

acceleration terms are not only kinematically related to the un~ 

• • * 

known input crank velocity 9 and acceleration 9, but also depend 
upon the elastic deformations, velocities and accelerations of 
deformation. Elastic deformations, velocities and the accelera- 
tions of deformation depend upon the inertia forces which in turn 
depend upon the acceleration terms. Since the acceleration terms 
and "the inertia forces are iater-dependent, this set of equations 
con be best solved iteratively with the set of equations (2.26), 
(2.35), (2.37), (2.38) aid (2.39). 

2.2.1 Iterative Schemes 

As stated earlier the equations of analysis alcngwith 
related equations cannot be solved directly and some iterative 
scheme is to bo adopted for their solution. The procedure adopted 
for thoir solution is shown on "the flow chart (Appendix A). The 
steps am as follows s 

gjjjp ^ j iho ''rigid” displacement, velocity and acceleration 

terms appearing cn the right hand side of "the 
equations (2.1 1) and (2,12) are expressed in terms 
• of velocity and acceleration of the crank. These 

quantities are then substituted in equations of 
analysis* 3J® start with, the links ©re considered 
t© be rigid. <■ 
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SEEP 2 : Equations (2,1l) and (2.12) can now be solved using 

some numerical technique (e.g. fourth, order Bunge - 
Kutta Method) . Since the above equations constitute 
a boundary - value problem, following procedure is 
adopted for their solution, 

(i) Following values are assumed in Hie first step* 

(a) Initial guess values for slider displacement 

! 

X, its derivative X with respect to inde- 
pendent parameter © and initial velocity of 

# *2 

the main drive 9 (or for Y = 9 ). 

(b) Applied torque constants a^ end ag* 

(ii) Using Bunge - Kutta method values of X, X md Y 
are computed for each 2° of input - crank rotation. 

(iii) The initial values of X^ X and Y at 9 = 0 and 
their values at 0 = 2/1" 316 compared. In the 
steady state they should be equal, for this de- 
notes that the mechanism maintains m average 
speed in operation. If X (9 = 0) ^ X (9 = 2J\ ) 
the initial guess forX at 9 = 0 is corrected. 
Similarly if x' (9 = o) ^ X* (9 = ZJ\) the ini- 

I 

tial guess value of X is corrected. If 

Y (9 = 0) ^ Y (9 = 2 7T ) then the term of 

*2 

the torque expression (T = a^ ~ Og 9 ) is sui- 


tably changed. 
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(iv) The average computed value of Y is compared 
with S L? where SI is the desired average 
angular velocity of crank. The average value of 
Y is obtained as 31 arithmetic mega of the 180 
computed values of Y based on equal step in ©. 

If average value, of Y is not close to -A the 
initial guess value of Y is modified. 

(v) If the requirements of steps (iii) aid (iv) are 

not satisfied, with the modified values of X, 

X* , Y and a* at 0 = 0 the equations (2.1l) aid 
7 1 

m 

(2.12) are again solved. The iterative scheme 
is continued till requirements of steps (iii) 
and (iv) are met. 

STEP 3 : Having solved equations (2.1l) and (2. 12) the initial 

position of the mechanism (© = o) is set. 

STJEP 4 : Solution of step ( 2 ) yields slider displacement, 

velocities and accelerations of the slider aid "the 
crank. Knowing accelerations the system forces are 
evaluated using expressions (2.30). 

STEP 5 : Elastic deformations may, then, be computed using 

expressions (2.26). 

STEP 6 1 Using oppressions (2.38 and 2.39) velocities and 

accelerations of deformation are computed. 
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SEEP 7 


STEP 8 


STEP 9 


STEP 10 


STEP 11 


s The values of acceleration terns are modified using 
results of step (6). 

* Values of new and old accelerations are compared. If 
the values do not match to the desired accuracy the 
entire procedure from step (4) and onwards is repeated. 
I'm ally the modified values of inertia forces and acce- 
lerations are obtained. 

s All the velocity terns are modified using results of 
step (6). Displacement of coupler end C is also mo di— 
fied using results of step (5). 

* The process from step (4) to step (9) is repeated for 
all mechanism positions between 0 and 2 JJ with a step 
increment of 2°. 

\ 

: If the procedure from step (2) to step (lo) is iterated 
once only, step (2) is 35am undertaken with modified 
displacements, velocities and accelerations and the 
procedure is repeated. 

: Hew and old values of slider displacement X are com- 
pared. If the values match to the desired accuracy 
the process is terminated. Otherwise the entire pro- 
cedure from step (2) to step (ll) is iterated. 


STEP 12 



CHAPTER III 


HESUIES MD DISCUSSIONS 


To illustrate the application of the mathematical model 
and solution technique developed in Chapter II three numerical 
examples have been worked out. All the xhree examples were taken 
from the published literature so that the results obtained may be 
compared with the existing ones. In the first example problem 
the mechanism considered by Davidson 5 J was studied. In other 
two example problems the mechanisms synthesized by liniecki jj26 
were worked out. 


5.1 Example 1 : 

The mechanism studied in this example probloa is sche- 
matically identical with that exa m i n ed by Davidson J~5 j • In 
Davidson f s model the radius of crank, length of connecting rod 
and the amount of offset are taken to be 2.68 in., 9.38 in. and 
0.0 in. respectively. The driving crank rotates at an average 
speed of 200 rev/min aid the slider weight is taken as 4 lb. wt». 
The stiffness of the spring is taken to be 42.5 lb/in,. The 
masses and elasticity of the crank and the connecting rod are 
neglected in his model. 'Whereas in the present dissertation lii^y 
are taken into account. The crank is assumed to have m area of 
cross - section of 1.0 sq. in. and moment of inertia about its 
axis of rotation as 1.0 lb-in-sec 2 . The coupler weighs 2.0 lb. 
and has an area of cross - section as 0.75 in. 


Its moment 



Angular velocity of crank 0 < rad/sec) 



Slider displacement X ( in.) 
Rigid displacement of coupler-end 



Fig. 6 Variation of slider displacement a 

displacement of coupler-end C with 
to crank angle 



Angular velocity of coupler cj> (rad /sec) 
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Fig. 7 Variation of angular velocity of coupler 
with respect to crank angle 


Velocity of slider X (in./sec) 

Linear velocity of CG of coupler'ln x-direction Xq( in./sec) 
Linear velocity of C G- of coupler m y-direction Yq (m./sec) 
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'Aflgtttfef acceleration of crank 8 {rad /sec) 
Angular acceleration of coupler ip (rad / sec^} 


■4 0-0 


Cran k angle 0 (Deg ) 


800 




V /" 


Fig, 9 Variation of angular acceleratfor 

crank and coupler with res^ej 
crank angle / 
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d 


4-0 



Elastic deformation at C-G of coupler in x-direction X„ (in x 10 ) 

Og 3 

Elastic deformation at CG. of coupler in y-direction y r ( in x 10 ) 


o 

o 


2-0 




-4*0 r 

Fig ]3 Variation of elastic deformations at centre o 
gravity of coupler in x and y-direc tions wit I 
respect to crank angle 


Inertia torque at the crank P| (Lb-m.) 
Inertia torque at C.G.of coupler P 0 (Lb-in.) 



Inertia force of coupler end C P/ ( Lb.) 


in 

CL 


i2-ei~ 


A p c 


u 

o au 


O © 

d d 


a> a> 
c c 


jr \ f \ 36i 


Crank angle 0 (Deg ) 


Fig. 15 Variation of Inertia forces at coupler end 
and at centre of gravity of coupler in x and 
directions with respect to crank angle 
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'iatior i of (X-X c ) with respect to crank 
jle 0 


DISPLACEMENT In 


Ad 


o 



i 


INPUT ROTATION 0 Deg 

Fig. 17 Motion characteristics for solution 
mechanism of Davidson's model 


ACCELERATION in^sec 


Deference between actual and kinematic velocities of slider ( m/se 



Fig. 18 Variation of difference between actual and 

Rigid velocities of slider with respect to crank 
angle (Examples 2and 3) 
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of inertia about the transverse axis passing through its centre of 
gravity is taken as 0.155 lb-in-sec. 2 . She modulus of elasticity of 
the material of crank and coupler is taken as 30 x I© 6 lb/in 2 . The 
constant 82 of the motor torque is taken as 1.8. ! 26l 

3*1*1 Results and Discussion For Example Problem 1 

With the above input data KEDA of the mechanism was carried 
out. The results obtained are represented graphic ally as shown in 
Figures 5 to 16. 

Figure 5 represents the fluctuation of angular velocity 

• 

© of the crank. The desired average speed of the crank is 21.0 
rad/sec while the actual average speed of -the crank is allowed to 
fall in the region represented in between the lines 1 and 2 (FIG-.5), 
permitting an error of 0.5$. 

Figure 6 represents the variation of slider displacement 
X and "rigid" displacement Z^ of coupler - end C with respect to 
craak rotation ©. In the absence of spring, and with rigid links, 
curve 2 would have represented the variation of slider displacement 
X. Therefore, the vertical ordinate between curves 1 and 2 at sny 
value of © represents the combined effect of inclusion of the 
spring and inherent elasticity of links of the mechanism on slider 
displacement Z at that position of the crank. It may be noted 

that the characteristics of the mechanism are calculated taking 

« 

into account the variation of crank speed © due to the changing 
inertia of the astern with crank rotation. 



51 


Figure 7 represents variation of angular velocity of the 

« 

coupler 0, TAfcile Figure 8 represents the variation of slider velo- 

# 

city X and x and y - directional velocities of the centre of gra- 
vity of the coupler with respect to 0. As nay be expected, the 
slider velocity is not zero at crank angles 0° s 180° and 360°. 

Figure 9 represents the variation of angular acceleration 

of crank aid coupler with respect to 0. Eefering to Figure 5} 

• * 

since the craik - speed decreases from. 0 = 0° to © = 48°, © is 

• • 

negative in this range. Between 0 = 48° and 0 = 130°, 0 is posi- 

tive as 6 continuously increases in this raage* Similarly between 

*# 

© = 130° and © = 275°, © is negative and from © = 275° to 0 = 340° 

* # • • 

0 is positive. From 340° to 360° of crank rotation, © is again 

negative. 


Variations of linear accelerations of centre of gravity 
of the coupler in x aid y-directions and slider acceleration with 
© axe shown in Figure 10. 


The variation of elastic deformation X at coupler end 

e 

C with 0 is represented in Figure 11, while Figure 12 end 13 re- 
present how the remaining element deflections vaiy with 0. The 

elastic defoliation X (FIG. 11 ) when added to 1 rigid’ displace- 
_ e 

ment X_ (FIG. 6) will yield the actual displacement X of the 
coupler end 0. The variation of the length of spring (X - _X^) 
with respect to © is shown in Figure 16, while the variations of 
inertia forces with © are plotted in Figure 14 and 15» 


| l V i 

6 * ^ * 1 !•** 1 * i 


Ar< 


25631 
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Befering to Figures 9, 10, 11, 14, 15 aid 16 it may be 

noticed that there is a ehsnge in the nature of the curves of 

• • *• 

©, X, X Q , P , P 4 and (X - X ) in between © = 160° aid © = 260°. 
e 

This may be attributed to the nature of variation of X^ aid 

e 

hence (X - X„) with respect to ©. The clastic deformation X 
o u e 

varies between 160° aid 260° of crank rotation in such a way that 
(X - X„) is affected in this range of 0 as shown in Figure 16. 

V 

The variation of (X - X Q ) directly affects inertia force P 4 

• * 

(equation 2.30) and slider acceleration X (equation 2.l), which 

*• 

in tum affect the input cr ack acceleration 0 and inertia torque 

The variations of slider displacaaaat, velocity and 
acceleration with input angle © as obtained by Davidson jjTj for a 
schematically identical mechanism, are shown in Figure 17. Compa- 
ring his results with those obtained in the present analysis it 
was found that as may be expected, they do not tally. lbs diffe- 
rence in the results of slider acceleration is much more pronounced 
in compansion to the differences in displacement and velocity 
results. As per Davidson's results the slider moves 4.0 in. as 
the crank rotates from 180° to 300° while the present analysis 
gives -this distance as 3.99 in. Davidson predicts the slider velo- 
city at © = 180° as 0.0 in./sec. , at © = 240° as 42.3 in/sec and 
at © = 300° as 66.0 in/sec. The corresponding values obtained by 

2,269, 39.352 aid 67.8 in/sec respectively. 


the present analysis are 
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Since Davidson did not consider the effect of inherent 
elasticity s the fluctuation in driving link speed and the mass of 

crank and coupler (viiidh have been accounted in the present analy- 

. \ ** 

sis; the present results for X are much different from one obtained 

•* 

by him. He obtains values of Z at 180° , 250° aid 300° of crank 

P 

rotation as 1900, 0.0 and - 900 in/sec approximately as against 

2 

964.4, 740.0 and - 300*0 in/sec obtained in the present analysis, 
3.2 Examples 2 and 3 

The mechanisms studied in examples 2 and 3 are schemar* 
tically identical to those synthesized by Idniecki j 26~j . In 
Liniecki's first model the radius of crank, length of coupler and 
the amount of offset are taken to be 5.889 in., 19.164 in. and 
0.688 in. while in the second model these values are taken as 
5*933 in., 26.961 in. and 2.59 in. respectively. In both the models 
following values have been adopted by Idniecki : 

-1 2 

Mass of the slider = 0.52 ibf in sec 

Moment of inertia 2 

of crank about its = 24.0 lbf in sec 

axis of rotation 

Driving Torque 

Constant =1.8 

Average speed of 

crank = 26.2 rad/sec 

Since Liniecki has considered rigidly attached slider, 
the stiffness of the spring in the present case is taken to be 
very high (52000.0 lbf/in.), Idniecki has not considered mass of 
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tile coupler separately (he considers equivalent moment of inertia 
of rotating masses and equivalent reciprocating mass in the mechar- 
nism) and has assumed "rigid” links. The following additional 
information is assumed for both the models examined here* 

Area of cross-section of crank = 2.0 in? 

Area of cross-section of coupler = 1*5 in? 

Modulus of Elasticity E = 30 * IQ 6 lbf/in? 

Moment of inertia of coupler 

about the transverse axis = 5.1 lbf in. sec? 

passing through its centre of 

gravity. 

The mass of the coupler is assumed to be half that of the slider. 

3.2.1 Results and Discussion of Example Problems 2 and 3 

The KEDA of the mechanisms of examples 2 and 3 was per- 
formed. The difference between the actual velocity of slider and 
the desired velocity fj26j was obtained for every 20° of crank 
rotation, Einally this quantity was plotted in Eigure 18 against 
crank rotation 0. 

Liniecki synthesizes the first mechanism by making the 
error between the actual and the desired velocities of the slider 
zero at 60° 3 80° and 120° of crank rotation. The second mechanism 
is synthesized by making this error zero at 60°, 160° and 240° of 
crank rotation. But the results of KEDA of both the mechanisms 
predict that the velocities at these precision points will not bo 
.zero (Pig4 18) if elasticity is also taken into account. This 
error at precision points varies from 3 % to 10$, 



CHAPTER IV 


CONCLUSIONS AND SCOPE FOR FUHEHER WORK 

The import aace of KEDA of mechanisms has long been 
re cognize a by the researchers. As it is clear from the example 
problems worked out, this analysis is essential for the mecha- 
nisms where accuracy is main objective, since the elastic 
deformations may cause appreciable error in mechanism charac- 
teristics. 

However, in the present dissertation the effects of 
clearances and tolerances, danping characteristics and friction 
in joints have not been taken into account. An analysis con- 
sidering all these effects together would be an ideal one and 
will predict the actual characteristics of the mechanism much 
more closely. Efforts on future may be directed to formulate 
a mathematical model of mechanisms considering all these effects 
and to device suitable solution techniques. 
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RFAD1 lu t ( ABC ( J) , J=*, y-r) 

100 FORMAT ( &F ID . 5 ) 

ABC(I9)=3000U^OU.-i 
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ABC ( x 3 ) =5 2,! . 

ABC( 14) =1 000. C 
12 CONTINUE 
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CALL ANALYS 

10 CONTINUE 
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DI MENS I ON DPHDOT { ivl ) » DXGjOT !bb) > C'YODuT { 1 £u ) » DC ETA ( lo * 3 t 'J^CC 
1181) , DACCYG ( iGl 3 tXGDs.1T ( !•- 3 t YGDOT ( 3 - 3 » XCFK 13*) 

CALL FLUNCOl.’ 3 

D A T — A j t^A TA M / I i \ 


AO ( 


T=PAI/90.0 
tPS5=5. 0 
EPS6-1.'* 
ZX(1)=ABC{1 3 
ZX(2)=ABC (2 ) 
ZX (3)=AB3 (3 3 
A1=ABC (4) 

R= ABC ( 5 3 
AL=ABC (6) 
A=ABC { T 3 
A2=ABC < 8 3 



01 =ASC ( 9 ) 

GI =ABC ( 1 .. ) 
om:=abcc:j 
GM 2 = ABC ( i ' ) 

S = 4 BC (lb) 

WN=ARC( 3 A) 

A 0 =ABC ( i 5 ) 

AG-ABCC.6) 

Ii=A 6 C ( :7 ) 

I 2 =ABC ( IS ) 
c = A3C (19) 

0 M£GA=ABC (BO) 

£P$ 1 =ABC( 21 ) 

2 PS 2 =ABC( 2 ? ) 

EPS 3 =A 3 C ( 22 ) 

£PS / -=ABC( >’■* ) 

IRc=, 

MM = ^ 

J .7 MM=*!M+i 
MK = 1 
NK = i 

DO i 1 = 1,2 
ZZ( 1 ,I)=ZX( I) 

NN=i 80 
N J = l*. 

N= iSi 
X= 0.0 
1111 K =1 

XCFF= 0 . : 

DO U' K J=x,,j 
XD ( K J ) = 0. 0 
KK=KJ+i 
Xl=X+T/ 2 . ’* 

X 2 =X+T 

CALL FUMCf X, ZX » CZX , K, XCFF ) 

AL PHA ( KJ ) =DZX ( 2 ) / 2 . j 
XDOUBL { KJ )*UZX{.-»> 

Y ( KJ ) =ZX( 2 ) 

ACCX(KJ)*ZX( 3 )*aLPH\{KJ)+Za<. }*,;c ’LBLtKJ) 
DO 20 1 = 1,3 

DZY 1 U)=DZXU ) 

23 ZY{I)=ZX( I)+. , 5 *T*DZX( I) 

CALL FUMC(X*,ZY,DZX,K».iCFF) 

DO 25 1 = 1,3 

0ZY2 ( I } =D ZX ( I ) 

25 ZYU)=ZXU)+ . 5 *T* 0 ZX( 1 ) 

CALL FUNIC( XI , ZY , QZ X , K, XC FF ) 

, DO 30 1 * 1,3 

.*, DZY 3 U)=DZXU> 

ZY(I)=ZX( I) +T*OZX(I ) 



0 1 = A 6 C { 9 ) 

G I = A B C ( 1 w ) 
tDM 1 = ABC ( 1 1 } 

0M2 = ABC ( i. 2 ) 

S=ABC { 13 ) 

WN=A8C( 34 ) 

AQ=A3C< 15 ) 

AG=ABC ( 16 ) 

U=ABC( 17 ) 

12 =ABC ( IB ) 
c= ABC( 19) 

OMEGA* ABC { 20 ) 

£PSI=A6C( 21 ) 

EPS2=ABC( 22 ) 

EPS3=ABC { 23 ) 

£PS4=ABC( 24 ) 

I Rc= . 

MM = y 

17 MM=MM+1 
MK=i 
NK = 1 

00 1 1=1,3 

1 ZZ(1,I )=ZX£ I ) 

NN = i80 
N J=l0 
N=18i 
X = 0.0 
Llli K= 1 

XCFF = 

DO i»' KJ = i., M 
XO { K J ) =0. 0 
KK=KJ+i 
Xl=X+T/2. J 
X2=X+T 

CALL FUMC ( X , ZX , C ZX , K t XC F F ) 

ALPHA (KJ)=DZX( 21/2.0 
XDOUBUKJ ) = U Z X { 3 ) 

Y ( KJ ) =ZX( 2) 

ACCX ( K J ) = ZX (3 ) #AL PH A { K J ) + ZX { 2 ) *X DOCBL { K J ) 
DO 2 < 1=1,3 

DZY1 ( I ) =DZX ( I ) 

2D ZY(I )=ZX( I)+0.5*T*DZX(I> 

CALL FUMC(Xi,ZY,DZX,K,XCFF) 

DO 25 1=1,3 

DZ Y2 ( I ) =D ZX { I ) 

25 ZY(1 )=ZX( I) +A.5*T*0ZX(1 ) 

CALL FUMC{Xi,ZY,DZX,K,XCFF) 

DO 30 1=1,3 

DZY3( I )=DZX ( I ) 

3D ZY ( I ) =ZX( I) +T*DZX(I ) 



CALL FUNC( X2, ZY, 11 ;,K,/.C“-F5 
DO 3 5 I=;,i 

DZY4( I )=DZX { I ) 

ZX ( I ) =ZX< I } +T*< DZY_ ( I ) + .. . *DZY . S 15 - t 2 . *_,ZY3 ( I 5 + DZYM 15 5/6. 
35 ZZ (KK, I ) = ZX ( I 5 
X= X+T 

10 CUNT I Nile 
b34 CONTINUE 

HHx = ZZ(l,l)-ZZ{.8„, .5 
HHZ = ZZ{1,25-ZZ{13:, :_) 

HH3=ZZ { It 3 5 - ZZ { 18*., 3) 

Y A V = i"' ,%j 

DO 797 I=_, _0 

797 YAV=Y4V+ZZ( I,Z) 

Y A V = Y A V / 1 SO . 0 
AV=OMEGA**2 
HH4=YAV-AV 

IF (ABS( HH1 5 .Lx. EPS! 5 70 TO ;50 
MK =2 

ZZ(i,l)=ZZ( 


250 

IF { A6S ( HH2) .Lb . r PS6 5 GO 
MK = 2 

A1=A1+5.0*HH2/3.C 

TO 

21 'Z' £- 

252 

IF t AS S ( HH3) .Lc.EPS3 JOG 
MK = 2 

ZZ(l,3)=ZZ(_,j ) -C .5 *HH2 

T 1/ 

25-t 

254 

IF IABSI HH4) .L2.SPS5 5 GO 
MK = 2 

ZZ C1,2)»ZZ( 1,2)-HH4 

T 

257 

25 7 

I F ( MK. NE. 2 ) GO TC ^5 
DO 795 1 = 1,3 
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ZX ( I ) =ZZ ( 1,1 5 

IF INK. EQ. 3) GO TC 335 




GO TO 17 



25 5 

CONTINUE 




PRINT99»MM 

99 FORMAT ( 5DX, *TGTAL I TF.RA7I UNS TAKEN**, Ti/SOX, 22 ( 1H-) ) 
PRINT750, Ai 

75 U FORMAT(20X, *A1=*,F15.3) 

PR INTI 00, (ZZ(K,1) ,K = I, 131, t ) ) 

100 FORMAT ( 50 X, *SL I DER DI S PLACEMENT* //6 { IUX, Fi 5. S, 5 (5X , F15 . b ) / ) 5 
PRINT101, (ZZ ( K, 2), K=l, 181,10) 

131 FORMAT ( 50 X, * VALUE OF Y*//6U..,X,F 5. 3 , 5 ( 5X , F i5 .8 ) / ) 5 
PR I NT102, ( Z Z ( K , 3 ) » K = l , 1 Ci , 1 ) ) 

132 FORMAT ( 50X, *VALUE OF XPR IM*//6 (a OX , Fx 5 . 3, 5 ( 5X , F 15. 8 ) / ) ) 

DO 200 1=1,181 

23a XDOT ( I )=ZZ( 1,3) *SQRT ( ZZ { I ,2) ) 

PRINT103, (XDOT(I), 1=1, 181,13) 

4.03 FORMAT ( 5DX, *VALUE OF XD0T*//6( IOX, Fi5. 3, 5 ( 5X, F15.3 ) / ) ) 
PRINT256 
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256 F0RMAT(5.*JX,*CJWH3* r -• : 5 ; - * ) 

833 CONTINUE 
1 =' 

X = ' . 

P3=R**2/( 4. C*uL ) 

P4=A**2/( 2. (.UAL ) 

P7=R*A/AL 
Pli=R/AL 
3 1 = 1 + ., 

I M I TA { I ) = 0 
11 = 1+1 
1 J=I +2 
XX =2. 0*X 

ANGVL'L { I ) =SyRl { Y C I ) ) 

P5=WN**2/ ( Wli** 2 -Y { I ) ) 

P6=WN**2/ . *Y { I } ) 

P 8 =AL*( 1. 0- { A/ AL + x*SI J( 5 / AL 1 )#*-. >. =5 
P9=PG**2 

Pi ■= ( A + R*SI N( X ) ) * (R*C03 ( X ) )** •' 

P12=R*SI N ( X )+ ■ . 5*Pi ( . 5*R*$IN ( XX)+A*CUS ( X) ) 
P13=R*C0S (X ) +0 • 5 # Pi * * ( R *C US { X X )-A*Ei j ( X ) ) 

Z = ARSIN(A/AL + R*SI’4( X)/AL) 

P14=E *AQ 
P15=E*I1 
P16 = £ ! XAG 

P17=SIN(X)+ TAN { Z) *C,JS(X ) 

P 1 8 = C OS ( X ) - T A ; J { Z ) *S I m ( X ) 

Pi9=C0S (X )-.5*TANCZ )*SIM X) 

P20=S IN { X ) + .5*TAN ( Z ) *C0 Si U 
P2i = TAN(Z )/COS(Z) 

P22*1.9/C0S(Z>**2 

P23=P13**2 

P24=P17**2 

88 l=C.5*AL**2*Pi7/Pi 5 

BB21=-SIN(X )*Pi 8 /(P i.4*CUS i l) ) 

B822=AL*Al.*C0S(X)*P17/( :-.*Pi 5*C0S( Z ) ) 

BB23=P21/P16 

8B2=8B21+8B22+BB23 

8B31=AL*P19*P1S/Pi*+.5*AL*P?2/P;6 

BB32=AL*AL*AL*P17*P20/{^.0*Pi5) 

BB3=-BB31-BB32 

BB4=AL*P23/PI4+AL*P22/P.16+AL**2*P?''>/{3.*P15) 

BB51=AL*PI8*SIN(X)/(2.*Px4) 

8B52=AL*AL*AL*PI7*CUS(X>/<6.*Pi5 ) 

6B5=8B52-BB5i 

XC{I)=R*C0S(X)+P8 

ALP 1=AL/ ( E*AQ ) 

ALP2=AL**3/ (3.0*E*I i> 

ALP3=AL**2/ (2-0*E*I 1) 

ALP4=AL/( E*I1) 
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*COS(Z ! ) 


•' L P 2 * ■ 
1 + ALP2*. 


5**2+AL P E 
S^^o+AL^? 
5*B?+ALPc 


ALP5=AL*P22/( -♦AG) 

ALP6=SI fi{ Z)*P22/{1:*4G) 

ALP7=-AL*P2 "/ ' 2 .0*t. *A0 3 

AL P8 = AL / U 2 . 0* ri* I 2 ) +7A ; , ( Z 3 ** 2/ { - L* 

ALP9=-0.5*ALPo 

ALP^-'-AL* (P?2+CGS(Z 3**1 )/ (<-. *:*A0 ) 
ALP1I=AL**3*SIN{ 2. *2 3/(91,. *=*I1)- 
ALPI 2=AL*SI f-i( Z )**_/ ( *.0*l *A0)+AL** . 
B1=-SIN (X )/{ AL*CCS( Z3 ) 

82=-P 19 
B3=Pi8 

B4 =- 0 • 5 *S I N ( X 3 
B5=COS ( X ) /{ *L*C0S(Z ! ) 

B6=-P20 
B7 = P i. 7 

68=0. 5*C3$( X) 

8B 6=ALP 3* 35 

887=ALPl*Bl**.i+AL?2*i'5**2 + ALP c 
BBS=ALPI*B1 *B2 + ALP2*:5*B6+AL ; >9 
869=ALPl*Bl*B3+ALP2* r 5*B7+ALPc 
BB 10=ALP1 *B 1*64+ ALP ..*B3*~ 3 
BB 1 1 = AL P3 *R 6 

BB12 = ALP1*31*32 + ALP2*B5*6 6+'4LP9 
BB13 = ALPl*Bl**2+ALP-,:*B6**-+ALPi.O 
BB 14= ALP1 *32*83+ ALP 2*B6*B'+ALP7 
BBI5 = ALPl*B2*B4+ALP 2*86*B-‘'+ALP.i- 
BB 16 = AL P3 *3 8 

BB17 = ALPl*Bl*84+ALPk*65*'3"> 

BBi8 = ALPl*B2*B4+ALPi:*B6*Hij+ALPlL 
BB 19 = AL PI *3 3*54 + ALP 2*87*3 0 
BB 20 = At PI *64**2 + ALP 2*58 **L + A LPi 2 
IF ( IRE.GE .1 )GU TO 962 
XCF<I )=XC (I ) 

CONTINUE 

XFOOT ( I ) =-ANGVLL ( I ) * ( R*P5 *SI N ( X ) +P 
ACCXGU ) = -P 12*ALPHA { I )-P12*Y(I 3 
ACCYG ( I ) = ;'•. 5*R* ( COS ( X ) * ALPHA (I >- SI 
BETA( I )=(R*CUS(X)*ALPHA(I ) + (Pl ./P9- 
Q1 ( I ) =OI*AL PHA ( I ) 

Q2 ( I ) =GI*BETA( I ) 

Q3 ( I )=0M2*ACCXG( I ) 

Q4(I )=S*( ZZ (I ,1 )-XCF( I ) 3 
Q5 ( I ) =0M2 *ACCYG ( I ) 

XCDEL ( I ) = BB I*QI ( I )+8B2*QZ ( I ) +863*0 
XCFU ) =XC ( I ) +XCDEL( 1 ) 

CONTINUE 

IMITA( I ) = IM IT At I )+I 

PHYOEL( I ) =BB6*Q1 ( I ) +8B7*Q2 ( I ) +B8 8* 

XGDEH I ) = BB11*QH I ) +BB 12*02 ( I 3+^6* 

YGDEL( I ) = B316*Q1 (I ) +BB17*Q2( I ) + «- B , 


) + LL**a*SI ; ( z )**? / ( . *-■*!/) 

-,.L*Sri{ i. *z )/( j. * 2* '■ G ) 

. *c jS< Z 3**2/ (4_.0*2 *Ia) 


' *P5*C r JS { X } + 2 .*P3*P6*S I ,M( XX 3 3 
X ) *Y ( I 3 ) 

-R*S IN { X ) ) *Y{ I ) )/Pd 


A I 3 +BB4*04 ( I ) + BB5*w5( I ) 


W3( I 3 +B89*Q4 { I) +BBl4*Q5 { 1 ) 
*Q3( I 3 +BB14*Q4( I )+BB15*Q5U ) 
■*Q3( I ) +5 r *’. 9* I J + BB20*Q5( 1 ) 
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IF ( I .EQ.M )3ti It: 5 5' 

HH i = { PHY3 EL (II )-PHY;. . L ( l) ) / T 
HH2= ( XGOE L( il J-XGDjUI ) }/7 
HH3= ( YGDEK I I ) -YGL‘1 L { I ))/T 
I F ( I.FQ.NM I J = :. 

HH*t= ( PH YD EL ( I J )-2 .0*PH Y J j:. ( I I } +? PY jL ( I) ) / 

HH 5= ( XGDEL( I J )-2 .(:•* \ ^ D L { i I ) + vr, tL( I) ) /T**: 

HH6= ( YGDt L( I J ) . *YGD:. L ( 1 1 ) +YuH £L ( I ) 5 / T * * Z 

DBcT A ( I ) =HH4*Y ( I ) +HH;.*£LPHA ( I ) 

DACCXC{ I) =HII5*Y{ I ) + i IH2* ,L { I ) 

DACCYG( I ) =HH6*Y { I ) + HH."'.* -\LPH4 ( I S 
GO TO 552 
551 CONTINUE 

D8 ETA ( I ) = D3ETM ]. ) 

DACCXG ( I ) =0 ACCXG ( 1 ) 

DACCYG( I)=DACCYG(i) 

55 2 CONTINUE 

6ETA{ I ) 5DBETA ( I ) + (R*COS (/) *ALPH* (I )+( P A 0/?9-?.*S I’i { X )) *Y ( I ) ) / ?.:. 

ACCXG ( I ) = -P 12*ALPHA ( I )-P„0*Y( I ) -r L L C X G { I 5 

ACCYG ( I )="DACCYG{ I)+ / . 5*8* ( CDS ( X } *Z^PHA{ I ) - S I j(X)*Y< I ) ) 

Q1 ( I ) =01* AL PHA { I ) 

Q2( I ) =GI*BETA ( I ) 

Q3(I)=0M2*ACCTG( I) 

Q4(I)=S*(ZZ(I,1)-XCF(I ) ) 

Q5(I )=0M2*ACCYG( I) 

XC DEL ( I ) = BB i*Q 1 ( I ) + NT 2 * 01 ( I ) + PB3 ( I ) +BB4*C4{ I ) + B3i*Q5 ( I) 

XCF1 ( I )=XC( I )+XCC! L (I } 

HHK= ( XCF1 { I l-XCF ( I ) 5/XCF( I ) 

XCF { I )=XCF1 (I ) 

I F { ABS ( HHK ) .LT . SP$4 ) GO TO 554- 
GO TO 555 
554 CONTINUE 

IF(I.EQ.M)GO TO 556 
DPHDOT ( I ) =HH1*ANGVEL ( I ) 

OXGDOT ( I) =HH2*ANGVEL ( I ) 

DYGDOT ( I ) =HH3*ANGVr L ( I ) 

GO TO 557 

556 CONTINUE 

DPHDOT ( I ) =DPHDOT ( 1 ) 

DXGDOT (I)=DXGD0T(1) 

DYGDOT ( I ) =DYGDOT ( 1 ) 

557 CONTINUE 

PHYDOTC I) ^DPHDOTt I)+Pli*COS(X )*ANG-/EL( I J/COS ( Z) 

YGDOT ( I ) =-DY GDUT ( I )+ ; . 5*R*C0S ( X ) * ANGVl L ( I ) 

XGDOT ( I )=HOXGDUT ( I )-AMGVEL { I) *<R*$I , (X)+0. 5*P 11* (0.5*R*SIN( XX) +A* 
ICOS(X))) 

X= X+T 

IF(I.LT.N) GO TO 3 

PRINT 638, (IMITAU), 1 = 1,181,10) 

688 FORMAT! 10 X, *NQ OF ITERATIONS BiiFCR- CONV-." RG .-\G2 */lX , 201 6 ) 
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20 i 
202 

203 

204 

205 

206 
207 

207 0 
212 

213 

214 

215 
215 2 

216 
217 
219 
22u 
221 
222 

223 

224 

225 

226 
22 7 
228 


PRINT 201, CXDD'JBL ( I ) j) ) 

FORMAT! 50 X, * VALUE OF vVii i f -.L *// *■ { j ;• , f ~ 
PRI.MT2C2, (ANGVuLt I), I = ) 

FORMAT { /5 .,X » *\f{ LU1 JF A JG/OL " *//w . /.s r 

PR I NT 203, ( A CC X ( I ) , I =;. , „ 

FORMAT (50X, *SL! DER -XC*//6( JX,F 
PR I NT2'.“ 4, ( XF7CJT { I ) , 1 = 1, .a. ,1 } 

FORMAT ( /5 IX , JF XFj_T *//:{ 

PR I NT 20 5, (A LPHA { I ) , 1 = j. , ,. y) 

FORMAT ( /50X , *ANG ACC*//o { . OX , F ' 5 
PRINT2J6, ( 8 FTA { I ) , I = ; , , _ j 

FORMAT ( /5 >X , *VALU E JP R -'T •- *//6{'‘ 

PR I NT 20 7, (ACCXGtl .,.;)) 

FORMAT ( /5 GX , *y ALU c OF iCC "G* / /r. ( : 

PRINT2J7D, ( ACCYG { I) , 1 = ^, . ) 


» 5 ( r. X , F „ 3 . 5 / 5 

*/ / 3 { „ X , = 3 F i . 1 

. ' , 3 ( 5 v , F .. 5 . 3 ) / ) ) 

X , F „ 5 . - , c l 5 X » F 5 . 2 ) / J ) 
,-:(px,f /) ) 

,F_i . 3( 5X, F «. 5 • ? )/) ) 

, r" f . ,'( s <, F .:.-)/) ) 


FORMAT ( /5 f X ,* aCCYG *//_,(■ X , F . i . , 5 { ~ X , F_ 5 . )/) 5 

PR I NT 2 12, (3i ( I ) , 1=1 , 0.,. » 1 >,: ) 

FORMAT (/5 OX, *VALU OF L‘..*//6( „;.X,F: s.O, - ( 5 < , F 15 . c }/) ) 
PRINT213, (U2(I),I=j ,5 0,.,, ) 

FORMAT ( /5 -»X,*VAttH: OF 01*/ /St: CM , F„ 5 . 2 , 5 { 5 X »F 15 . 3 ) / } ) 
PRINT214, (Q3 { I ) , 1 = 1 , 1 Ox , 1 0 ) 

FORMAT ( /5 OX, *V4Ll_ JF 0 2 * / / 6 ( 1 G X , F * 5 • 2 , 5 ( f X , F i 5 . i) / ) ) 
PRINT215, (Q4( I ) , 1=1, _ 51,. } 

FORMAT ( / 5 ;X , *VALUc OF 0'*//G! 1CX,F15. 2 , 5 ( 5 X , F 15 . a ) / ) ) 

PR I NT2152 » ( Q5 ( I ! , I = ! , 101, .. 0) 

FORMAT ( 50X, *VALUE UF CG*//6( 2 X, Flf . 2 , 5 { 5 V . , FI 5. .3 ) / ) ) 
PRINTS 16, (XCOMLt I), 1 = 1, 18.,,! ,) 

FORMAT ( /5 0X ,* VALUE UF LCDcL * / / 6 ( ± U X , FI 5 .t ,5 ( 5 X , F 1 5 . 

PRINT217, (XC ( I ) , 1 = 1, •- 81,10} 

FO RMAT ( /5 ■ X , *VALUF OF XC*/ ft ( CX , F „5 . 0 , 8 { 5X , F 15 . 3 } / ) ) 


PRINT219, { XCF ( I ) , 1= i , 13 * , . •) 

FORMAT t/5 OX, *F I HAL XCF*//6( t OX , FI 0. e , 5 ( 5X, F.5 . 8 ) / ) ) 
PRINT220, (PMYOELd ) ,1=1 ,1U, iO) 

FORMA T ( 5vX, *PHYD EL*//6 { 1 X, F i*> . ? , 5 { 5X , F 15 . ? ) / 5 ) 

PR I NT 221 , (XGDifL ( I ) , 1 = 1, *e. ) 


FORMAT ( 50 X, *XGDcL*//6{ 1 OX , FI 5.8, 5( 5X , FI 5. 3)/) ) 


PRINT222, ( YGUEL ( I), 1 = 1, 181, 10) 

F0RMAT(5DX,*YGD6L*//6( 1 ;X,F15.S, 5( 5X, Fi5-C) /) ) 
PRINT223, { PHYDOT ( I ) , 1 =1 , 1 31 , 1 

FORMAT { 50X, *PHYDCT*//6{ iOX,F15. 8 »5 (5> ,F15.8) / ) ) 


PRINT224, (XGDOTt I ), 1=1,131,10) 

FORMAT ( 53 X,*XGD0T*//6(1 X, FI 5. 6 , 5C 6X, F 1 5 . 8 ) / > ) 
PRINT225, (YGDOT ( I ), 1 = 1, lb 1,1” ) 

FORMAT (50X,*YGD0T*//6(10X,F15.B,5( \,Fi5.F)/) ) 
PRINT 226, (DPHDOTC I ) ,1=1,131 , x s) 

FORMAT ( 51 X, *DPHD0T*//6( ), j X , F 15. 8 ,5 { 5X , F 15 . 8 ) / ) ) 
PR I NT 227, ( DXGDOT (I), 1 = 1, 181,1 ) 

FORMAT ( 50 X, *0XGDCT*//6 ( 10X, F 15 . 8 ,5 ( 5X , F15 . 8 ) / ) ) 

PRINT228, (DYGOOT(I) , 1=1,181, 10) 

FORMAT t 50 X, *DYGDCT*//6 ( X, F 15. 8 ,5 ( 5X , F 15 .8 ) / ) ) 


3 

)/) ) 


6 )/)) 



PRINT229, ObcTA(I),I = : , .. , >) 

229 FORMAT ( 53X, *DE2TA*//r. ( , F. 5. X, F = .•’)/)> 

PR 1 NT 23** » ( DaCCXG { I } j I = '; , , ' ’ ) 

230 FORMAT150X, *DtCCX>//6?. r { y ,F:5. )/)) 

PRINT231, (DACCYSU) , 1 = .-., . ) 

231 FORMAT ( 50X, *DACCYG*//e- 1 . F .5. ■; ,C i K \ , F .3 . > / } ! 

DO 837 1 = 1 ,i.5L 

HK=( ZZ( I, 1) -XD{ I ) )/Ll { I , 1 ) 

I F { A B S ( HK ) . G _ . 0 . 1 ; ... ) GO T'j 

837 COOT I NUc 
PR I NT 33 9 

839 FORMAT { 53 X,*PR0BL;.-1 HAS CO.-.*/ - KG- ..*/ SOX, I’ - ( - H* ) ! 
RETURN 

838 DO 840 1 = 1, iS. 

840 XD { I ) = ZZ ! 1,1) 

IR t= I Rt+1 
PRINT961, IR c 

961 FORMAT ( /5 OX , *R1F IT 10=*, I •!/ 5 X , *1 i J // ) 

MM = u 

836 CONTINUE 
MM=MM+1 

K=2 
X=3 . 0 

DO 809 1=1,3 

809 ZX(I)=ZZU,I) 

DO 8 lu JK=i,N 
KK=JK+i 
KK J=KK+1 
XI =X+ 0 . 5*T 
X2=X+T 
XCFF=XCF( JK) 

XCFP= (XCDtL (KK)-XCD :..L( JK) )/T 
ZS ( 1 ) =PHYDOT ( JK ) 

ZS ( 2 ) =XGDOT { JK ) 

ZS(3)=YGD0T (JK) 

ZS ( 4) =8ET A{ JK) 

ZS ( 5 ) =ACCXG { JK ) 

ZS ( 6 ) =ACC YG ( JK ) 

CALL F UN C ! X,ZX,CZX»K,XCFF) 

ALPHA! JK) =DZX(2)/2.0 
XDOUBL t JK ) = DZX ( 3 ) 

Y { JK)=ZX( 2) 

ACCX! JK)=ZX (3)*ALPHA{ JK)+ZX( 2)*XCGUBL( JK) 

DO 820 1=1,3 

DZ Y1 ( I ) =DZX ! I ) 

820 ZY( I )=ZX( I) +0.5*T*DZX!I ) 

XCFF=0.5*(XCF{JK)+XCF(KK) ) 

XCFP= r 3« 5* (XCFP+ ( XCDcL ! KKJ ) -XCDtL !KK ) } /T ) 

ZS ( 1 ) = .* . 5*( PHY DOT ! JK ) +PHYDOT ( KK ) 5 
ZS (2) =0.5*< XGDOT! JK )+XGDOT IKK) ) 
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ZS (3 )=•.. .5*( YGDUT ( JK ) +YG 337 { KK ! ) 

ZS(4)=-;.E3*(KFT4(JK) :fT { K K ) ) 

ZS ( 3 ) =0 .5 -'*( ACCXu ( JK ) + •'> C C / ; { K ?> ) ) 

ZS (6) = J.5*( ACCYGt JK)+;c:Y;'.(K<) ) 

CALL FUNC ( X; , Z Y , CZ X , K , / C FF ) 

DO 8 25 I = _ , j, 

DZ Y2 ( I ) =DZX { I ) 

825 ZY ( I ) = ZX( I ) + .. . 5*T**JZX ( I ) 

XCFF=u.5*(XCF( JK)+XCF(KK] ! 

CALL FUMC ( Xi ,ZY, DZ,\,X , -C P F) 

DO 830 1=1,3 

DZ Y3 ( I ) =D ZX ( I ) 

83 L ZY C I )=ZX( I)+T*CZX(I ) 

XCFF=XCF{ KK ) 

XC FP= ( XC3EL ( KKJ )-XC J 1 L UK ) ) / T 
ZSU)=PHYDOT{KK) 

ZS ( 2 ) =XGDOT ( KK } 

ZS { 3 ) =YGDOT ( KK ) 

ZS (4) =BcT A{ KK ) 

ZS(5)=ACCXG (KK) 

ZS (6 ) =ACCYS (KK ) 

CALL FUN C ( X2 , Z Y , OZ < , K , x C F F ) 

DO 835 1=1,3 

DZY4( I ) =DZX ( I ) 

ZX ( I )=ZX( I) +T*(DZY1( I ) + 2 . *D ZYi ( I ) +2 . Y" ( I ) + D Z Y 4 ( I ) ) /6.U 

835 ZZ(KK,I)=ZX(I) 

X=X+T 

810 CONTINUE 
MK=3 
MK=3 

GO TO 834 
END 
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*IBFTC FUUC 


SUBROUTINE FUNG ( X t Z j 3Z j K , ,-.C FF 5 

* if * # * * * * * $ * * * * jj, £ ^ ^ . V ^ .,. ^ .j, ^ . lf . # ^ ^ ^ 


THIS SUBROUTINE EV-iL'iATrS 
APPLYING RUNG'Z-L! *TT , 
COMMON/TA J/ZS 
CQMMON/SMITA/XCFP 
COMMON/KHARt/R, AL»A,;'I, il 
DIMENSION ZS( 6) 

DIMENSION Z(i),CZ('.) 

XX = 2. *X 


THE F J r.C . I . . v; L..i REQ'jIT 
■ETHGL FE : THE 3 ~ L ; i 7 I ' J F 


rJ’G » •' ” 2 j E f iv ! i r A 1 » A E 


PI = { A+R#S IN ( X 5 ) / AL 


P2=A/AL 


P3 = R/AL 


P4 = P1 

P5 = i.i'-P4 

P6=SI N ( X) **2 

P7 =COS { X ) *#2 

P3=CnS(X) **3 

P9=Ul+GI*P3**2*P7/P5 

P10 = 0M2 5!c R**2*{ P6+0. l**P'(* P4+ 0.3*51 ( KY ) *P1 ) 
SSi=P9+PID 

Pli =GI*P3**2* ( 2. >J*P3*P8*P.-S IN ( '< X ) P 5 )/P5 

Pi2=0M2*R**2*( SIN ( XX )+ . 5*P? *PS*P^-0. 25*S I M ( XX ) *P4> 

Pi3=QM2*R**2*(C0S(X<)*P.. + u„5*P2*SI , ( X X ) *COS(X 5 3 

SS2=P11+P12+P13 

SS3=-R*C3S(X)-AL*P5**. . 5 

SS4=R*SIN (X )+R*CCS< X)*P. 

SS 5 = WN*#2 *( R*GOS (X)+AL*P5**0.S 3 
H5u=UM2*R**2/4. ,• 

H51=H5D*C0S (X)**2 


H52=H5 )*S IN (XX) 
SS 1=SS1 +H51 
SS2 = SS2-H 5? 


DZ ( 1 ) =Z ( 3 ) 

P14=A1-S* (Z (1)+SS3)*(Z(3)+SS* )-Z (2 3 * { A2+0 .5*SS2 ) 
DZ(2)=2.0*( P14-CJ M#Z(3)*(SS5-Z(1 )*Wi**d 3 } /SSL 
DZ (3) = (SS5-Z(1 )*KN**2-G.5*Z( 3}*0Z( .) 3/ZC2) 

IF (K.EQ.l 3 RETURN 
SS4=SS4-XCFP 
SS3=-XCFF 
SS5=WN**2*XCFF 

P14=A1-S*(Z (l)+SS3)*(Z(i)+SSH)-Z (2 )*A2 

Pi 5=G I *ZS ( i )*ZS{4)+DM2*(ZS(2)*ZS (5 3 +ZS ( 3 3 *ZS ( 6 ) 3 

P16=-P15/ SQRT ( Z ( 2 3 3 

P17=0M1*Z (3 ) « ( Z ( 1 3 *WN**2-SS5 3 

DZ (23=2.3* ( Pi 4+P 16+ P 173 /OI 

DZC3)=<SS5-Z<1 )*WN**2-0.5*Z(3)*DZU) )/Z(2) 

RETURN 

END 


W H I L ; 
:ALYSi 
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